Excluded- volume effects in the diffusion of hard spheres 
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Excluded-volume effects can play an important role in determining transport properties in diffu- 
sion of particles. Here, the diffusion of finite-sized hard-core interacting particles in two or three 
dimensions is considered systematically using the method of matched asymptotic expansions. The 
result is a nonlinear diffusion equation for the one-particle distribution function, with excluded- 
volume effects enhancing the overall collective diffusion rate. An expression for the effective (collec- 
tive) diffusion coefficient is obtained. Stochastic simulations of the full particle system are shown to 
compare well with the solution of this equation for two examples. 
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I. INTRODUCTION 

Recently there has been an increasing interest in un- 
derstanding the transport of particles with size-exclusion 
Size exclusion is important in many biological pro- 
cesses, including diffusion through ion channels [1, Q and 
in chemotaxis and can have a significant impact on 
the thermodynamics and kinetics of biological processes 
such as association reactions at membranes Finite- 
size effects are also important when considering the com- 
bustion of powders [6|, collective behavior (e.g. animal 
flocks or traffic movement) 0,13 and granular gases @. 

Excluded-volume or stcric interactions arise from the 
mutual impenetrability of finite-size particles (see Fig. 
[T]). For one-dimensional configurations, such as chan- 
nels, the single-file diffusion of hard-core particles can be 
solved exactly by mapping it to the classical diffusion 
of point-particles [l^, HH. This has recently been ex- 
tended to heterogeneous particles and anomalous parti- 
cles However, the situation in higher dimensions 
is more challenging. 
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FIG. 1. (Color online) Excluded [red (dark grey) and black] 
and available [blue (light grey)] area in a solution of black 
particles for the placement of an additional test particle, (a) 
The area available with point particles is the whole domain, 
(b) The area available (to the center of the test particle) with 
finite-size particles is reduced. Modified from Minton |14 ]. 



It is well-known that for finite-size particles the 
effective diffusion coefficient becomes concentration- 



dependent. In fact we have to distinguish between two 
alternative notions of diffusion coefficient: the collec- 
tive diffusion coefficient, which describes the evolution 
of the total concentration, and the self-diffusion coeffi- 
cient, which describes the evolution of a single tagged 
particle • Here we concentrate on collective diffusion, 
hereafter simply referred to as diffusion. 

Batchclor [l^ models Brownian diffusion of particles 
with hydrodynamic interaction using generalized Ein- 
stein relations to find a concentration dependent correc- 
tion to the collective diffusion coefficient. Felderhof [T3l 
considers the same problem through an analysis of the 
Fokker-Planck equation, and includes both excluded vol- 
ume and hydrodynamic effects. His analysis is based on 
the thermodynamic limit (in which the number of parti- 
cles N and the system volume V tend to infinity, with 
the concentration N/V fixed), and is valid only for a 
small perturbation from the equilibrium concentration. 
Similarly, the self-diffusion coefficient to first-order in a 
constant concentration is obtained from the generalized 
Smoluchowski equation in [H, [T^ . 

Muramatsu and Minton [l^ use a simple model to cal- 
culate the diffusion coefficient of hard spheres by esti- 
mating the probability that the target volume for a step 
in a random walk is free of any macromolecules. Other 
authors model excluded volume phenomenologically by 
introducing a particle pressure, resulting in an equation 
of state in which the c omp ressibility is reduced as the 
concentration increases [20| . 

Another popular approach is to consider lattice mod- 
els, in which a particle can only move to a site if it is 
presently unoccupied. Such an approach has been used 
to model diffusion of multiple species with size exclu- 
sion effects [2f|, [l^ or to model the effect of crowding on 
diffusion-limited reaction [2^. 

The preceding approaches are all either phcnomcno- 
logical in nature, restricted to small perturbations from 
a uniform concentration, or based on the thermodynamic 
limit in which the number of particles tends to infinity. 
Here we consider a finite number of finite-sized particles 
diffusing in a box of fixed size. We perform an asymp- 
totic analysis of the associated Fokker-Planck equation 
in the limit that the volume fraction of particles is small. 
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Our analysis is systematic, using the method of matched 
asymptotic expansions, but is not appropriate for con- 
centrations close to the jamming limit. 

II. DIFFUSION WITH FINITE-SIZE EFFECTS 

In order to focus on steric effects, we suppose that there 
are no electrostatic or hydrodynamic interaction forces 
between particles. We work in d dimensions, where d is 
either 2 or 3. Thus our starting point is a system of N 
identical hard core diffusing and interacting spheres (or 
disks) , each with constant diffusion coefficient Dq and di- 
ameter , in a bounded domain ft in M.'^ of typical diam- 
eter L. By nondimensionalizing length with L and time 
with L'^/Dq, the size of the domain and the diffusion coef- 
ficient may be normalized to unity, while the diameter of 
the particles becomes e = K/L. We assume that the par- 
ticles occupy a small volume fraction, so that Ne'^ <C 1. 
We denote the centers of the particles by Xi(i) g at 
time t > 0, where 1 < i < N Each center evolves 

according to the stochastic differential equation (SDE) 

dX, = \/2 dB, + f, dt, l<i<N, (I) 

where the are N independent d-dimensional stan- 
dard Brownian motions and is the external force on 
the i th particle. In general this force may include both 
inter-particle and external interactions, such as electro- 
magnetic, friction, convection and potential forces, in 
which case depends on the positions of all the particles 
X = (Xi, . . . , X^r). While soft-core steric effects can also 
be built into fi, hard-core collisions can be more easily 
expressed as reflective boundary conditions on the "col- 
lision surfaces" r = ||X, — Xj|| = e, with 1 < i < j < N. 
Since we are focusing on hard-core particle interactions, 
we restrict ourselves to other external forces of the form 

F(X) = [f(Xi),...,f(X^)], (2) 

where f : ^ R'' acts identically on all N particles. 
We suppose that the initial positions Xi(0) are also ran- 
dom, and that they arc independent and identically dis- 
tributed. 

Let P{x, t) be the joint probability density function of 
the N particles. Then, by the ltd formula, P{x, t) evolves 
according to the linear Fokker-Planck partial differential 
equation (PDE) 

dP ^ ^ 

-^=Vs-[VsP-F{x)P] in (3a) 

where V g and • respectively stand for the gradient 
and divergence operators with respect to the A^-particle 
position vector x = (xi,...,XAr) e Vt'^ . Note that be- 
cause of steric effects, ([3a| is not defined in $7^ but in its 
"hollow form" f]f = n^\B^, where = {x € VL^ -.^i ^ 
j such that ||xi — Xjll < e} is the set of all illegal con- 
figurations (with at least one overlap). On the collision 



surfaces we have the reflecting boundary condition 

[Vj P - F{x) P]-n = on dfl^ , (3b) 

where n G S'^^~^ denotes the unit outward normal. 
Since the initial positions of the particles are indepen- 
dent and identically distributed, the initial distribution 
function Po{x) is invariant to permutations of the par- 
ticle labels. The form of ([3]) then means that P itself 
is invariant to permutations of the particle labels for all 
time. 

Although linear, the PDE model ^ is very high- 
dimensional, and it is impractical to solve it directly. 
Since all the particles are identical, we are interested 
mainly in the marginal distribution function of the first 
particle, given by p{xi,t) = J P{x,t) dx2...dxAr. We 
aim to reduce the high-dimensional PDE for P to a low- 
dimensional PDE for p through a systematic asymptotic 
expansion as e — )■ 0. 



A. Point particles 

In the particular case of point-particles (e = 0) the 
model reduction is straightforward. In this case the N 
particles are independent and the domain is D,^ = VL^ 
(no holes) , which implies that the internal boundary con- 
ditions in (j3b[) vanish. Therefore P{x,t) = YliLiPi^i^^)^ 
and 

dv 

^(xi,<) = -[Vx.p-nxOp] in n, (4a) 
0= [VxiP-f(xi)p]-ni on dn, (4b) 

where fii is the outward unit normal to dft. Note that 
since the particles are indistinguishable each satisfies the 
same diffusion equation and boundary condition, so that 
P is a product of N identical I-particle distribution func- 
tions p. If the particles were not identically distributed 
initially then we would need a different distribution func- 
tion for each one; although these would all satisfy the 
same diffusion equation they would have different initial 
conditions. This point will be important when we go on 
to consider finite-sized particles. 



B. Finite-size particles 

When e > 0, the internal boundary conditions in (j3b[) 
mean the particles are no longer independent. When we 
integrate ([5a|) over X2, . . ., xjv and apply the divergence 
theorem we end up with surface integrals over the colli- 
sion surfaces, on which P must be evaluated. However, 
when the particle volume fraction is small, the volume 
in occupied by configurations in which three or more 
particles are close is small [0{e'^'^N'^)] compared to those 
in which two particles alone are in proximity [0{e'^N)]. 
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Thus the dominant contribution to these "colhsion in- 
tegrals" corresponds to two-particle collisions. We illus- 
trate our approach for = 2; since two- particle collisions 
dominate the extension to arbitrary N is straightforward. 
A similar approach is used in [l^ [l^ ■ 

For two particles at positions Xi and X2 , Eq. ((3a|) reads 



dP 

'at 



(xi,X2,t) = Vxi • [VxiP-f(xi)P] 



Vx, •[Vx,P-f(x2)P], (5a) 



for (xi,X2) G fi^, and the boundary condition (|3b[) reads 

[Vx,P - f (xi)P] • fii + [Vx,P - f (X2)P] • n2 = 0, (5b) 

on Xi G dil and ||xi— X2II = e. Here = ni/||ni||, where 
rii is the component of the normal vector ft corresponding 
to the i—th particle, n = (ni, n2). We note that fii = 
on X2 G dQ, and that fii = —112 on ||xi — X2II = e. 

We denote by ri(xi) the region available to particle 2 
when particle 1 is at Xi, namely r2(xi) = \ Be(xi). 
Note that when the distance between Xi and dfl is less 
than e the volume |0(xi)| increases. This creates a 
boundary layer of width e around dfl where there ex- 
ists a wall-particle-particle interaction (three-body inter- 
action). Since the dimensions of the container are much 
larger than the particle diameter these interactions are 
higher-order and we may safely ignore them and take 
|0(xi)| constant [2^. Integrating Eq. ([5a|) over il(xi) 
yields 



dp 



xi,t)= Vxi • [VxiP-f(xi)p] 



/ [f(xi)P-2VxiP- Vx,P]-n2d52 (6) 
[Vx,P-f(x2)P]-n2d52. 



anuaB,(xi) 

The first integral in ^ comes from switching the order 
of integration with respect to X2 and differentiation with 
respect to Xi using the transport theorem; the second 
comes from using the divergence theorem on the deriva- 
tives in X2. Using (j5bp and rearranging we find 

^(xi,t)= Vx, •[Vx,p-f(xi)p] 

+ / {-2Vx,P + P[f(xi)-f(x2)]}-n2d52. (7) 

Because the pairwise particle interaction is localized near 
the collision surface 5Pe(xi) we can determine it using 
the method of matched asymptotic expansions [2^ . 



C. Matched asymptotic expansions of the density P 

We suppose that when two particles are far apart 
(||xi— X2II 3> e) they are independent, whereas when 



they are close to each other (||xi — X2II ^ e) they are 
correlated. We designate these two regions of the con- 
figuration space f2g the outer region and inner region, 
respectively. 

In the inner region, we set Xi = xi and X2 = xi + ex 
and define P(xi,x, t) = P(xi,X2,<). With this rescaling 
^ becomes 

BP 

e2_(xi,i,t) = 2V|P-62Vx, • [f(xi)P] +e2v|^P 
+ eVx • { [f (xi) - f (ii + ei)] P} - 2eVi, • ViP, (8a) 
with 

2i • ViP = ex • {Vi,P + [f(ii + ei) - f(ii)] P}, (8b) 

on ||x|| = 1. As noted above, we can neglect the bound- 
ary layer and hence assume that Xi is not close to dft; 
the region in which the particles are close to each other 
and the boundary is even smaller, and will affect only the 
higher-order terms. In addition to (|8b[) the inner solution 
must match with the outer solution as x — > cx). In the 
outer region, by independence, 

P(xi,X2,t) = q{xi,t)q{x2,t), 

for some function q{x,t). Note that the invariance of P 
with respect to a switch of particle labels means that in 
the outer region both particles have the same distribution 
function q (as happened in the point particle case). The 
normalization condition on P gives q{x.i,t) = p(xi,t) -|- 
0{e'^). Expanding this outer solution in inner variables 
gives 



P(xi , X2 , t) = (?(xi , t)q{ii + ex) 

^ g^(xi,t) +eg(xi)x- Vii(?(xi) 



(8c) 



Expanding P in powers of e, and solving ((8a|) . (j8bp with 
the matching condition (|8cp we find that the solution in 
the inner region is simply 



P(xi,x,i) ^ q (xi,t) + eq(xi)x • Vi^^q{xi) 



(9) 



Using this solution to evaluate the integrals in ^ we 
find, to 0{e'^), 



dp 
di 



(xi,t) = Vx, • [Vxi(p + «rfeV) -f(xi)p] , (10) 



where a2 — 7r/2 and = 2tt/3. The extension from two 
particles to particles is straightforward up to 0{e'^), 
since at this order only pairwise interactions need to be 
considered. Particle 1 has {N — 1) inner regions, one 
with each of the remaining particles. A similar procedure 
shows that the marginal distribution function satisfies 

^ (xi , i) = Vx, • { Vx, [p + ad{N- l)e V] - f (xi ) p} . 

(11) 

Equation (jlip describes the probability distribution func- 
tion for finding the first particle at position xi at time 
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t. Since the system is invariant to permutations of the 
particle labels, the marginal distribution function of any 
other particle is the same. Thus the probability distri- 
bution function for finding any particle at position xi at 
time t is simply Np. 



these corrections is that due to interactions between three 
(or more) particles. Because our asymptotic expansion is 
systematic, these correction terms could in principle be 
calculated. 



D. Interpretation 



III. COMPARISON WITH THE FULL 
PARTICLE SYSTEM 



We see from (fTTj) that steric interactions lead to a 
concentration-dependent diffusion coefficient, with the 
additional term proportional to the excluded volume. 
Equation ([TT|) is consistent with that derived by Felder- 
hof [l7| . but extends it to situations in which p is not 
close to uniform. We emphasize that (fTTj) is valid for any 
N. However, for large TV such that N — 1 N wc can 
introduce the volume concentration c — TTNe'^p/2d and 
rewrite ^ as [13 

dc 

-(xi, - Vx, • {D{c) V.,c - f (xi) c}, (12) 

where D{c) is the concentration-dependent collective dif- 
fusion coefficient, given by 



D{c) = 1 -|-4(d- l)c. 



(13) 



Note that the collective diffusion coefficient D{c) is in- 
creased relative to point particles. This is in contrast 
to the self-diffusion coefficient (which may be related to 
the mean squared displacement of a particular particle) 
which is reduced relative to point particles [l5| . This ap- 
parent contradiction may be understood as follows: the 
diffusion of any particular particle is impeded by its colli- 
sions with other particles. However, these collisions bias 
the random walk towards areas of low particle density, 
so that the overall spread of all particles is faster. To 
analyze the self-diffusion coefficient in the current frame- 
work we would need to label a particular particle, rather 
than treating all particles as identical. 

Whilst the self-diffusion coefficient can be thought of 
a diffusion coefficient intrinsically attached to each parti- 
cle, the collective diffusion coefficient relates the diffusive 
flux to the concentration gradient of all particles [2^ : the 
distribution function p is the probability of finding any 
particle at a given position, rather than the probability 
of finding a particular particle there. Thus the collective 
diffusion coefficient is not associated with an individual 
(tagged) particle or even a representative particle. This 
also means that it cannot easily be related to the mean- 
squared displacement of particles. This distinction has 
important consequences when upscaling from individual 
to collective behavior. 

In ([TT|) we have only included the leading-order non- 
linear term due to steric effects. There will be correction 
terms of 0{e'''~^^N) due to higher-order terms in the two- 
particle inner solution as well as new inner regions 
where three particles [©(e^'^iV^)], or two particles and the 
boundary [0{e'^^-^N)], are close. The most important of 



In order to assess the validity of ([TT]) we compare its 
solution p(xi,t) (obtained by a simple finite difference 
method) with Monte Carlo (MC) simulations of the 2N- 
coupled SDE ([T]) in two dimensions. The particle-particle 
(and particle- wall) overlaps are treated as in [2^ . To test 
the importance of steric interactions, we also compare 
with the corresponding solutions with e = 0. 
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FIG. 2. (Color online) Marginal distribution function p(xi,t) 
at time t — 0.05 with normally distributed initial data and 
N — 400. (a) Solution p(xi,t) of Q for point particles (e = 
0). (b) Histogram for e = 0. (c) Solution p(xi,t) of for 
finite-sized particles (e = 0.01). (d) Histogram for e — 0.01. 
Histograms computed from 10* realizations of H]) with At = 
10~^. All four plots have the same color bar. 



In Fig. [2] we show the results of a time-dependent sim- 
ulation with i = 0, n ^ [-i, e = 0.01, N = 400, 
for which the initial distribution is a Gaussian of zero 
mean and standard deviation 0.09 (normalized so that 
its integral over 51 is one); the figures correspond to time 
t = 0.05. The simulation time-step At is chosen such 
that almost no collisions are missed. The theoretical pre- 
dictions for both point and finite-size particles compare 
very well with their simulation counterparts, while steric 
effects are clearly appreciable even though the volume 
fraction of particles is only 0.0314. However, note that 
while the average concentration is low, the local concen- 



tration is considerably high at the origin: c = 0.617 at 
time t = and c = 0.0479 at time t = 0.05. The initial 
profile, in which particles arc concentrated in the center, 
spreads faster when steric effects are included [Fig. ^c)] 
than when they are not [Fig. HJa)], indicating that the 
overall collective diffusion is enhanced. 

When the force field f is the gradient of a potential, 
f(xi) = — VxiV^(xi), we may write PT|) as 

g = v., • (pu), (14) 

with u = Vxi [\np + 2adiN - l)e'^p + V{:xii)]. Equa- 
tion (fT4)) has an associated free energy [sOl 

F{p)= [ plnp + ad(iV-l)eVdxi + / F(xi)pdxi, 
Jn Jn 

where the first integral corresponds to the internal energy 
and the second integral is the potential energy. Note 
that excluded- volume effects increase the internal energy 
of the system. The stationary distribution, which we 
denote Ps (xi ) , is obtained by minimizing the free energy 
or by solving 

lnp,(xi) + 2ad{N - l)e'^p,(xi) + V{xi) = C, (15) 

with the constant C determined by the normalization 
condition on ps . For our second example we consider the 
volcano-shaped potential V(xi) = -4.77 e"^°°ll''ill' + 

^nl 111^ 

3.58 e""' "''^" and we compare the stationary distribu- 
tion Ps predicted by with simulations using the 
Metropolis-Hastings (MH) algorithm [3l|. Figure [3] 
shows the model and simulation results with N = 1000 
and ft and e as in Fig. [2] for both point and finite-size 
particles (with volume fraction 0.079 and volume con- 
centration c = 0.189 at the origin). In this case there is 
competition between the potential well and steric repul- 
sion: the particle density inside the well is reduced for 
finite-size particles. Again, the agreement between the 
model ([15]) and the stochastic simulations is excellent. 

IV. DISCUSSION 

We have derived systematically a nonlinear diffusion 
equation which describes steric interactions in the limit 
of small but finite particle volume fraction. Our method 
justifies for example the ansatz made in [3^ to account 
for the finite size of the cells in the Keller-Segel model 
and prevent aggregation, and unlike [itI |33| does not rely 
on a closure assumption. 

The equation we have derived is for the one-particle 
distribution function, which measures the probability of 
finding any particle at a given position; the particles we 
consider are identical and indistinguishable. This means 
that we are examining collective diffusion, and we find 
that this is enhanced by the finite size of the particles. 
We have not considered the self-diffusion of a particular 
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FIG. 3. (Color online) Stationary marginal distribution func- 
tion ps(xi) under the external potential V for point particles 
and finite-size particles, with A'^ = 1000. (a) Point particles, 
Ps (X . (b) Histogram for e = 0. (c) Finite-size particles Ps 
from (|15p (e = 0.01). (d) Histogram for e — 0.01. Histograms 
computed with 10^ steps of the MH algorithm. All four plots 
have the same color bar. 



(tagged) particle, which can be related to an individual 
particle's mean-square displacement. To analyze the self- 
diffusion coefficient in the current framework we would 
need to label a particular particle, rather than treating 
all particles as identical; we intend to do this in a future 
work in which we consider multiple particle populations. 

We note that for point particles in one dimension 
(where particles must move in single file and are not al- 
lowed to pass) HI] has observed density dependence in 
the self-diffusion (mean square displacement) of a tagged 
particle. Their interpretation is that the expansion of the 
whole system from dense to dilute environments quick- 
ens the self-diffusion of any tagged particle. This is a 
different effect to the one we have observed, since, as 
mentioned in the introduction, the collective diffusion of 
point particles in one dimension is linear, with a diffusion 
coefficient independent of the density. In two dimensions 
self-diffusion is less sensitive to the particle density, since, 
informally, there is much more space for particles to pass 
each other. 

The method we have developed was implemented here 
in its simplest setting (hard-core identical spherical par- 
ticles with an external potential) but it can be extended 
in many directions. Particles of different size (c/. [1^) or 
shape can easily be incorporated, while the hard-core in- 
teraction between particles can be replaced by any short- 
range soft-core interaction. 

On the other hand, incorporating long range effects 
such as chcmotaxis or electrostatic interactions is more 
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challenging; in such cases a system size expansion is likely 
to be needed in addition to a small particle expansion. 
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